﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SystemAnalizer.Common;

namespace SystemAnalizer.Models
{
    public class DiscreteMaker
    {
        private double m_CenterOfApproximation;
        private double H_step;
        private int m_DiscretesQuantity;
        private string m_Parameter = "t";
        private string m_Epression = "";

        public int DiscretesQuantity
        {
            get { return this.m_DiscretesQuantity; }
            set { this.m_DiscretesQuantity = value; }
        }

        public DiscreteMaker(double tnyu, double H, int K, string parameter)
        {
            setCenterOfApproximation(tnyu);
            setScaleCoefficient(H);
            setDiscretesQuantity(K);
        }


        //calculating discretes
        public double[] calculate(string expr)
        {
            setExpression(expr);
            int pos = 0;
            return parseExpression(m_Epression, ref pos);
        }

        //setting valid expression
        private void setExpression(string expr)
        {
            if (this.m_Parameter != "t")
                expr.Replace(m_Parameter, "t");

            expr = expr.Replace(" ", "");
            expr = expr.Replace("exp", "0$");
            expr = expr.Replace("sin", "0~");
            expr = expr.Replace("cos", "0?");
            expr = expr.Replace("tan", "0\"");
            expr = expr.Replace("ctan", "0\'");
            expr = expr += char.MaxValue;
            this.m_Epression = expr;
        }

        //set members
        private void setDiscretesQuantity(int K)
        {
            m_DiscretesQuantity = (K >= 1) ? K : 1;
        }

        private void setScaleCoefficient(double H)
        {
            this.H_step = (H >= 1) ? H : 1;
        }

        private void setCenterOfApproximation(double tnyu)
        {
            this.m_CenterOfApproximation = (tnyu >= 0) ? tnyu : 0;
        }

        //parse expression
        //parsing both operands of addition  or subtraction
        private double[] parseExpression(string str, ref int pos)
        {
            double[] result;

            result = parseMultiplyDivision(str, ref pos);
            while (str[pos] != char.MaxValue)
            {
                char op = str[pos];
                if (op != '+' && op != '-')
                    return result;
                ++pos;

                double[] right = parseMultiplyDivision(str, ref pos);
                if (result == null || right == null)
                    return null;
                if (result.Length > 0 && right.Length > 0)
                    if (op == '+')
                        result = DAddition(result, right);
                    else
                        result = DSubtract(result, right);
                else
                    result = null;
            }

            return result;
        }

        //parsing both operands of multiplying  or division
        private double[] parseMultiplyDivision(string str, ref int pos)
        {
            double[] result;
            result = parsePowExp(str, ref pos);
            while (str[pos] != char.MaxValue)
            {
                char op = str[pos];
                if (op != '*' && op != '/')
                    return result;
                ++pos;

                double[] right = parsePowExp(str, ref pos);
                if (result == null || right == null)
                    return null;

                if (result.Length > 0 && right.Length > 0)
                    if (op == '*')
                        result = DMultiply(result, right);
                    else
                        result = DDivision(result, right);
                else
                    result = null;
            }

            return result;
        }

        //parsing both operands of power or exponent
        private double[] parsePowExp(string str, ref int pos)
        {
            double[] result;
            result = parseSinCos(str, ref pos);
            while (str[pos] != char.MaxValue)
            {
                char op = str[pos];
                if (op != '^' && op != '$')
                    return result;
                ++pos;

                double[] right = parseSinCos(str, ref pos);
                if (result == null || right == null)
                    return null;
                if (result.Length > 0 && right.Length > 0)
                    if (op == '^')
                        result = DPow(result, right);

                    else
                        result = DExp(right);
                else
                    result = null;
            }
            return result;
        }

        //parsing sin or cos
        private double[] parseSinCos(string str, ref int pos)
        {
            double[] result;
            result = parseTanCtan(str, ref pos);
            while (str[pos] != char.MaxValue)
            {
                char op = str[pos];
                if (op != '~' && op != '?')
                    return result;
                ++pos;

                double[] right = parseTanCtan(str, ref pos);
                if (right == null)
                    return null;
                if (right.Length > 0)
                    if (op == '~')
                        result = DSin(right);
                    else
                        result = DCos(right);
                else
                    result = null;
            }
            return result;
        }

        //parsing tan or atan
        private double[] parseTanCtan(string str, ref int pos)
        {
            double[] result;
            result = parseCoefficient(str, ref pos);
            while (str[pos] != char.MaxValue)
            {
                char op = str[pos];
                if (op != '\"' && op != '\'')
                    return result;
                ++pos;

                double[] right = parseCoefficient(str, ref pos);
                if (result == null || right == null)
                    return null;

                if (result.Length > 0 && right.Length > 0)
                    if (op == '\"')
                        result = DTan(right);
                    else
                        result = DCTan(right);
                else
                    result = null;
            }
            return result;
        }

        //parsing coeficcient
        private double[] parseCoefficient(string str, ref int pos)
        {
            double[] result;
            bool negative = false;

            if (str[pos] == '-')
            {
                negative = true;
                ++pos;
            }

            if (str[pos] == '(')
            {
                ++pos;
                result = parseExpression(str, ref pos);
                if (str[pos] != ')')
                    result = null;
                ++pos;
            }
            else
            {
                string token = "";
                while (char.IsLetterOrDigit(str[pos]) || str[pos] == '.')
                {
                    token += str[pos];
                    ++pos;
                }

                if (token == "t")
                    result = DT_m(1);
                else
                {
                    try
                    {
                        double ans = System.Convert.ToDouble(token);
                        result = DMultiply(ans, DTed());
                    }
                    catch (Exception e)
                    {
                        result = null;
                    }
                }

                int i;
                if (negative)
                    if (result != null && result.Length > 0)
                        for (i = 0; i <= m_DiscretesQuantity; i++)
                            result[i] = -result[i];
            }
            return result;
        }

        // necessary elementar functions for discretes calculation
        // 1-s image
        private double[] DTed()
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            int i;
            Z[0] = 1;
            for (i = 1; i <= this.m_DiscretesQuantity - 1; i++)
                Z[i] = 0;
            return Z;
        }

        //  t^m image
        private double[] DT_m(int m)
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            double[] ted = DTed();
            double temp;
            int i, k;
            for (k = 0; k < m_DiscretesQuantity; k++)
            {
                Z[k] = 0;
                for (i = 0; i <= m; i++)
                {
                    if (k - i < 0)
                        temp = 0;
                    else
                        temp = ted[k - i];

                    Z[k] += MathHelper.binom(m, i) * Math.Pow(m_CenterOfApproximation, m - i) * Math.Pow(H_step, i) * temp;
                }
            }

            return Z;
        }

        // exp(X(t)) image
        private double[] DExp(double[] X)
        {
            double[] Y = new double[this.m_DiscretesQuantity];
            int k, i;
            Y[0] = Math.Exp(X[0]);
            for (k = 0; k < m_DiscretesQuantity - 1; k++)
            {
                Y[k + 1] = 0;
                for (i = 0; i <= k; i++)
                    Y[k + 1] += (i + 1) / (k + 1) * X[i + 1] * Y[k - i];
            }
            return Y;
        }

        //sin(x(t)) image
        private double[] DSin(double[] X)
        {
            double[] Y = new double[this.m_DiscretesQuantity];
            double[] Z = new double[this.m_DiscretesQuantity];
            int k, i;
            Y[0] = Math.Asin(X[0]);
            Z[0] = Math.Acos(X[0]);
            for (k = 0; k < this.m_DiscretesQuantity - 1; k++)
            {
                Y[k + 1] = 0;
                Z[k + 1] = 0;
                for (i = 0; i <= k; i++)
                {
                    Y[k + 1] += (i + 1) / (k + 1) * X[i + 1] * Z[k - i];
                    Z[k + 1] += -(i + 1) / (k + 1) * X[i + 1] * Y[k - i];
                }
            }

            return Y;
        }

        //cos(x(t)) image
        private double[] DCos(double[] X)
        {
            double[] Y = new double[this.m_DiscretesQuantity];
            double[] Z = new double[this.m_DiscretesQuantity];
            int k, i;
            Y[0] = Math.Asin(X[0]);
            Z[0] = Math.Acos(X[0]);
            for (k = 0; k < this.m_DiscretesQuantity - 1; k++)
            {
                Y[k + 1] = 0;
                Z[k + 1] = 0;
                for (i = 0; i <= k; i++)
                {
                    Y[k + 1] += (i + 1) / (k + 1) * X[i + 1] * Z[k - i];
                    Z[k + 1] += (i + 1) / (k + 1) * X[i + 1] * Y[k - i];
                }
            }
            return Z;
        }

        //tg(x(t)) image
        private double[] DTan(double[] X)
        {
            return DDivision(DSin(X), DCos(X));
        }

        //arctg(x(t)) image
        private double[] DCTan(double[] X)
        {
            return DDivision(DCos(X), DSin(X));
        }
        //ln(x(t)) image
        private double[] DLn(double[] X)
        {
            double[] Y = new double[this.m_DiscretesQuantity];
            double[] Z = new double[this.m_DiscretesQuantity];
            int k, i;
            if (X[0] == 0) throw new Exception();
            Y[0] = Math.Log(X[0]);
            Z = DDivision(DTed(), X);
            for (k = 0; k < this.m_DiscretesQuantity - 1; k++)
            {
                Y[k + 1] = 0;
                for (i = 0; i <= k; i++)
                    Y[k + 1] += (i + 1) / (k + 1) * X[i + 1] * Z[k - i];
            }

            return Y;
        }

        //D - options
        // multiply
        private double[] DMultiply(double[] X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            int i, k;

            for (k = 0; k < this.m_DiscretesQuantity; k++)
            {
                Z[k] = 0;
                for (i = 0; i <= k; i++)
                    Z[k] += X[k - i] * Y[i];
            }

            return Z;
        }

        private double[] DMultiply(double X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];

            for (int k = 0; k < this.m_DiscretesQuantity; k++)
                Z[k] = X * Y[k];

            return Z;
        }


        private double[] DDivision(double[] X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            double S;

            int i, k;

            if (Y[0] == 0)
                throw new DivideByZeroException();

            Z[0] = X[0] / Y[0];

            for (k = 1; k < this.m_DiscretesQuantity; k++)
            {
                S = 0;
                for (i = 1; i <= k; i++)
                    S += Z[k - i] * Y[i];
                Z[k] = (X[k] - S) / Y[0];
            }

            return Z;
        }

        //convert X(tnyu) to X(tnyu+1)
        //double[] DXtMyu(double[] X)
        //{
        //    int i;
        //    double[] Z = new double[this.m_DiscretesQuantity];

        //    Vector_double Xn_1(m_DiscretesQuantity, 1);
        //    Vector_double Xn(m_DiscretesQuantity, 1);

        //    Matrix_double P(m_DiscretesQuantity, m_DiscretesQuantity);

        //    for (i = 0; i < Xn.GetRow(); i++)
        //        Xn.SetElem(i, 0, X[i]);

        //    P.Zeros();

        //    for (i = 0; i < P.GetRow(); i++)
        //        P.SetElem(i, i+1, i+1);

        //    double newCenterOfApproximation = H;
        //    double tau;

        //    tau = (newCenterOfApproximation - m_CenterOfApproximation)/H;
        //    P = P*tau;
        //    P.ExpA();

        //    Xn_1 = P * Xn;

        //    for (i = 0; i < m_DiscretesQuantity; i++)
        //        Z[i] = Xn_1.GetElem(i, 0);

        //    return Z;
        //}

        //convert X(tnyu+1) to X(tnyu)
        //double *D::DXtNyu(const double *X)
        //{
        //    int i;
        //    double *Z = new double[m_DiscretesQuantity];

        //    Vector_double Xn_1(m_DiscretesQuantity, 1);
        //    Vector_double Xn(m_DiscretesQuantity, 1);

        //    Matrix_double P(m_DiscretesQuantity, m_DiscretesQuantity);

        //    for (i = 0; i < Xn_1.GetRow(); i++)
        //        Xn_1.SetElem(i, 0, X[i]);

        //    P.Zeros();

        //    for (i = 0; i < P.GetRow(); i++)
        //        P.SetElem(i, i+1, i+1);

        //    double newCenterOfApproximation = H;
        //    double tau;

        //    tau = (newCenterOfApproximation - m_CenterOfApproximation)/H;
        //    P = P*(-tau);
        //    P.ExpA();

        //    Xn = P * Xn_1;

        //    for (i = 0; i < m_DiscretesQuantity; i++)
        //        Z[i] = Xn.GetElem(i, 0);

        //    return Z;
        //}

        private double[] DAddition(double[] X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            for (int k = 0; k < this.m_DiscretesQuantity; k++)
                Z[k] = X[k] + Y[k];
            return Z;
        }

        private double[] DSubtract(double[] X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];
            for (int k = 0; k < this.m_DiscretesQuantity; k++)
                Z[k] = X[k] - Y[k];
            return Z;
        }

        // Y(t)^X(t) image
        private double[] DPow(double[] X, double[] Y)
        {
            double[] Z = new double[this.m_DiscretesQuantity];

            int k;
            bool isNumber = true;
            bool isParameterOnly = true;
            for (k = 1; k < this.m_DiscretesQuantity; k++)
                if (Y[k]!=0)
                    isNumber = false;

            double[] Ztmp = new double[this.m_DiscretesQuantity];
            Ztmp = DT_m(1);
            for (k = 0; k < m_DiscretesQuantity; k++)
            {
                if (X[k] != Ztmp[k])
                    isParameterOnly = false;
            }

            if (isParameterOnly && isNumber)
                Z = DT_m((int)Y[0]);
            else
                Z = DExp(DMultiply(Y, DLn(X)));

            return Z;
        }
    }
}
